home *** CD-ROM | disk | FTP | other *** search
/ Software 2000 / Software 2000 Volume 1 (Disc 1 of 2).iso / education / e077.dms / e077.adf / ScienceDemos / J2000 (.txt) < prev    next >
AmigaBASIC Source Code  |  1991-07-12  |  9KB  |  328 lines
  1. ' Program  "J2000"
  2.  
  3. ' copyright (C) 1986 by David Eagle
  4. ' 7952 W. Quarto Dr., Littleton, CO 80123, (303) 972-4020
  5. ' released into the public domain on March 20, 1986
  6.  
  7. ' conversion of stellar positions, proper motions, parallax and
  8. ' radial velocity from the standard epoch B1950 (FK4) to J2000 (FK5)
  9. ' reference: Astronomical Almanac, 1985, page X
  10.  
  11. OPTION BASE 1
  12. DEFDBL a-z
  13.  
  14. DIM a1(3),a2(3),upv(3),uvv(3),r1(3),r2(6),r3(6),v2(3),m(6,6)
  15.  
  16. SCREEN 1,640,200,3,2
  17. WINDOW 5,"Program J2000",(0,0)-(630,185),0,1
  18. PALETTE 4,0,0.8,0.2:' green
  19. PALETTE 5,1,1,0:' yellow
  20. PALETTE 6,0.8,0,0.93:' purple
  21. PALETTE 7,0.93,0.2,0:' red
  22.  
  23. pi=5.30795e-315
  24. pidiv2=pi/2
  25. deg.to.rad=pi/180
  26. rad.to.deg=180/pi
  27. xe=pi/12
  28. xer=12/pi
  29.  
  30. a1(1)=-5.19985e-315
  31. a1(2)=-5.18759e-315
  32. a1(3)=-5.18149e-315
  33. a2(1)=0.001244
  34. a2(2)=-0.001579
  35. a2(3)=-0.00066
  36.  
  37. ' define elements of transformation matrix
  38. m(1,1)=5.29981e-315
  39. m(1,2)=-5.26578e-315
  40. m(1,3)=-5.25963e-315
  41. m(1,4)=5.20278e-315
  42. m(1,5)=-5.16936e-315
  43. m(1,6)=-5.16294e-315
  44. m(2,1)=5.26578e-315
  45. m(2,2)=5.29981e-315
  46. m(2,3)=-5.22096e-315
  47. m(2,4)=5.16936e-315
  48. m(2,5)=5.20278e-315
  49. m(2,6)=-5.12435e-315
  50. m(3,1)=5.25963e-315
  51. m(3,2)=-5.22095e-315
  52. m(3,3)=5.29981e-315
  53. m(3,4)=5.16294e-315
  54. m(3,5)=-5.12434e-315
  55. m(3,6)=5.20278e-315
  56. m(4,1)=-0.000551
  57. m(4,2)=-0.238565
  58. m(4,3)=0.435739
  59. m(4,4)=5.29981e-315
  60. m(4,5)=-5.26578e-315
  61. m(4,6)=-5.25963e-315
  62. m(5,1)=0.238514
  63. m(5,2)=-0.002667
  64. m(5,3)=-0.008541
  65. m(5,4)=5.26578e-315
  66. m(5,5)=5.29981e-315
  67. m(5,6)=-5.22097e-315
  68. m(6,1)=-0.435623
  69. m(6,2)=0.012254
  70. m(6,3)=0.002117
  71. m(6,4)=5.25963e-315
  72. m(6,5)=-5.22095e-315
  73. m(6,6)=5.29981e-315
  74.  
  75. ' initialization
  76. CLS
  77. select%=1
  78.  
  79. CALL yesno.menu("Program introduction ?",intro%)
  80. IF intro%=1 THEN CALL introduction
  81.  
  82. ' main program loop
  83. WHILE select%=1
  84.  CLS
  85.  COLOR 3,0
  86.  PRINT "Right ascension with respect to 1950"
  87.  PRINT "< hours (0-23), minutes (0-59), seconds (0-59) >"
  88.  INPUT ras.hr.1950,ras.min.1950,ras.sec.1950
  89.  PRINT
  90.  COLOR 4,0
  91.  PRINT "Declination with respect to 1950"
  92.  PRINT "< degrees (-90 to +90), minutes (0-59), seconds (0-59) >"
  93.  INPUT dec.deg.1950,dec.min.1950,dec.sec.1950
  94.  PRINT
  95.  COLOR 5,0
  96.  PRINT "Right ascension proper motion with respect to 1950"
  97.  PRINT "< arc seconds per tropical century >"
  98.  INPUT ras.pmotion.1950
  99.  PRINT
  100.  PRINT "Declination proper motion with respect to 1950"
  101.  PRINT "< arc seconds per tropical century >"
  102.  INPUT dec.pmotion.1950
  103.  PRINT
  104.  COLOR 6,0
  105.  PRINT "Parallax with respect to 1950 < arc seconds >"
  106.  INPUT parallax.1950
  107.  PRINT
  108.  COLOR 7,0
  109.  PRINT "Radial velocity with respect to 1950 < kilometers per second >"
  110.  INPUT radialv.1950
  111.  
  112. ' convert 1950 right ascension and declination to radians
  113. ras.1950=xe*(ras.hr.1950+ras.min.1950/60+ras.sec.1950/3600)
  114. dec.1950=deg.to.rad*SGN(dec.deg.1950)*(ABS(dec.deg.1950)+(dec.min.1950*60+dec.sec.1950)/3600)
  115.  
  116. ' calculate unit position vector
  117. upv(1)=COS(ras.1950)*COS(dec.1950)
  118. upv(2)=SIN(ras.1950)*COS(dec.1950)
  119. upv(3)=SIN(dec.1950)
  120.  
  121. k1=21.095*radialv.1950*parallax.1950
  122.  
  123. ' calculate unit velocity vector
  124. uvv(1)=-ras.pmotion.1950*SIN(ras.1950)*COS(dec.1950)-m2(1)*COS(ras.1950)*SIN(dec.1950)+k1*upv(1)
  125. uvv(2)=ras.pmotion.1950*COS(ras.1950)*COS(dec.1950)-m2(1)*SIN(ras.1950)*SIN(dec.1950)+k1*upv(2)
  126. uvv(3)=dec.pmotion.1950*COS(dec.1950)+k1*upv(3)
  127.  
  128. ' E-term constants
  129. k2=a1(1)*upv(1)+a1(2)*upv(2)+a1(3)*upv(3)
  130. k3=a2(1)*upv(1)+a2(2)*upv(2)+a2(3)*upv(3)
  131.  
  132. FOR i%=1 TO 3
  133.  r1(i%)=upv(i%)-a1(i%)+k2*upv(i%)
  134.  r2(i%)=r1(i%)
  135.  v2(i%)=uvv(i%)-a2(i%)+k3*upv(i%)
  136.  r2(i%+3)=v2(i%)
  137. NEXT i%
  138.  
  139. ' matrix multiplication
  140. FOR i%=1 TO 6
  141.  a=0
  142.  FOR j%=1 TO 6
  143.   a=a+m(i%,j%)*r2(j%)
  144.  NEXT j%
  145.  r3(i%)=a
  146. NEXT i%
  147.  
  148. rm.2000=SQR(r3(1)^2+r3(2)^2+r3(3)^2)
  149.  
  150. ' compute declination wrt j2000
  151. CALL arcsin(r3(3)/rm.2000,b)
  152. CALL convert(b*rad.to.deg,dec.deg.2000,dec.min.2000,dec.sec.2000)
  153.  
  154. ' compute right ascension wrt j2000
  155. CALL atan3(r3(2),r3(1),c)
  156. CALL convert(c*xer,ras.hr.2000,ras.min.2000,ras.sec.2000)
  157.  
  158. k4=r3(1)^2+r3(2)^2
  159.  
  160. ' compute proper motions wrt j2000
  161. ras.pmotion.2000=(r3(1)*r3(5)-r3(2)*r3(4))/k4
  162. dec.pmotion.2000=(r3(6)*k4-r3(3)*(r3(1)*r3(4)+r3(2)*r3(5)))/(rm.2000^2*SQR(k4))
  163.  
  164. ' compute parallax and radial velocity wrt j2000
  165. IF parallax.1950=0 THEN radialv.2000=radialv.1950 :ELSE radialv.2000=(r3(1)*r3(4)+r3(2)*r3(5)+r3(3)*r3(6))/(21.095*parallax.1950*rm.2000) 
  166. parallax.2000=parallax.1950/rm.2000
  167.  
  168. CLS
  169. CALL display.data
  170. CALL yesno.menu("Another selection ?",select%)
  171. WEND
  172.  
  173. WINDOW CLOSE 5
  174. SCREEN CLOSE 1
  175. END
  176.  
  177. SUB convert(deg,hr,min,sec) STATIC
  178. ' convert degrees to hms or dms subroutine
  179. hr=FIX(deg)
  180. c=ABS(hr-deg)*60
  181. min=INT(c)
  182. d=(c-min)*60
  183. sec=INT(d+0.5)
  184. END SUB
  185.  
  186. SUB display.data STATIC
  187. ' display data subroutine
  188. SHARED ras.hr.1950,ras.min.1950,ras.sec.1950,dec.deg.1950,dec.min.1950,dec.sec.1950,ras.pmotion.1950,dec.pmotion.1950,parallax.1950,radialv.1950
  189. SHARED ras.hr.2000,ras.min.2000,ras.sec.2000,dec.deg.2000,dec.min.2000,dec.sec.2000,ras.pmotion.2000,dec.pmotion.2000,parallax.2000,radialv.2000
  190. WINDOW 4,"J2000 Data                   (press left mouse button to continue)",(0,0)-(630,185),0,1 
  191. CLS
  192.  PRINT
  193.  COLOR 1,0
  194.  PRINT TAB(25);"* 1950 <FK4> *"
  195.  a$=STR$(ras.hr.1950)+" hours"+STR$(ras.min.1950)+" minutes"+STR$(ras.sec.1950)+" seconds"
  196.  COLOR 2,0
  197.  PRINT TAB(5);"Right ascension";TAB(60-LEN(a$));a$
  198.  a$=STR$(dec.deg.1950)+" degrees"+STR$(dec.min.1950)+" minutes"+STR$(dec.sec.1950)+" seconds"
  199.  COLOR 3,0
  200.  PRINT TAB(5);"Declination";TAB(60-LEN(a$));a$
  201.  PRINT
  202.  COLOR 4,0
  203.  PRINT TAB(5);"Right ascension proper motion";
  204.  LOCATE ,48
  205.  PRINT USING "########.###";ras.pmotion.1950
  206.  COLOR 5,0
  207.  PRINT TAB(5);"Declination proper motion";
  208.  LOCATE ,48
  209.  PRINT USING "########.###";dec.pmotion.1950
  210.  PRINT
  211.  COLOR 6,0
  212.  PRINT TAB(5);"Parallax        (arc seconds)";
  213.  LOCATE ,48
  214.  PRINT USING "########.###";parallax.1950
  215.  COLOR 7,0
  216.  PRINT TAB(5);"Radial velocity (kilometers per second)";
  217.  LOCATE ,48
  218.  PRINT USING "########.###";radialv.1950
  219.  PRINT
  220.  COLOR 1,0
  221.  PRINT TAB(25);"* 2000 <FK5> *"
  222.  a$=STR$(ras.hr.2000)+" hour"+STR$(ras.min.2000)+" minutes"+STR$(ras.sec.2000)+" seconds"
  223.  COLOR 2,0
  224.  PRINT TAB(5);"Right ascension";TAB(60-LEN(a$));a$
  225.  a$=STR$(dec.deg.2000)+" degrees"+STR$(dec.min.2000)+" minutes"+STR$(dec.sec.2000)+" seconds"
  226.  COLOR 3,0
  227.  PRINT TAB(5);"Declination";TAB(60-LEN(a$));a$
  228.  PRINT
  229.  COLOR 4,0
  230.  PRINT TAB(5);"Right ascension proper motion";
  231.  LOCATE ,48
  232.  PRINT USING "########.###";ras.pmotion.2000
  233.  COLOR 5,0
  234.  PRINT TAB(5);"Declination proper motion";
  235.  LOCATE ,48
  236.  PRINT USING "########.###";dec.pmotion.2000
  237.  PRINT
  238.  COLOR 6,0
  239.  PRINT TAB(5);"Parallax        (arc seconds)";
  240.  LOCATE ,48
  241.  PRINT USING "########.###";parallax.2000
  242.  COLOR 7,0
  243.  PRINT TAB(5);"Radial velocity (kilometers per second)";
  244.  LOCATE ,48
  245.  PRINT USING "########.###";radialv.2000
  246.  WHILE MOUSE(0)=0:WEND
  247.  WHILE MOUSE(0)<>0:WEND
  248.  WINDOW CLOSE 4
  249. END SUB
  250.  
  251. SUB atan3(sinangle,cosangle,angle) STATIC
  252. ' four quadrant inverse tangent subroutine
  253. SHARED pidiv2
  254. IF ABS(sinangle)=0 THEN
  255.  angle=(1-SGN(cosangle))*pidiv2
  256. ELSE
  257.  angle=(2-SGN(sinangle))*pidiv2
  258. END IF
  259. IF ABS(cosangle)>0 THEN angle=angle+SGN(sinangle)*SGN(cosangle)*(ABS(ATN(sinangle/cosangle))-pidiv2)
  260. END SUB
  261.  
  262. SUB arcsin(sinangle,angle) STATIC
  263. ' inverse sine subroutine
  264. SHARED pidiv2
  265. IF ABS(sinangle)>=1 THEN
  266.  angle=SGN(sinangle))*pidiv2
  267. ELSE
  268.  angle=ATN(sinangle/SQR(1-sinangle^2))
  269. END IF
  270. END SUB
  271.  
  272. SUB yesno.menu(request$,response%) STATIC
  273. ' yes/no request subroutine
  274. WINDOW 4,request$,(0,0)-(215,45),0,1
  275. LOCATE 1,1:PRINT PTAB(25);"press left mouse"
  276. LOCATE 2,1:PRINT PTAB(25);"button to select"
  277. LINE (20,20)-(80,40),1,b
  278. LINE (140,20)-(190,40),1,b
  279. LOCATE 4,1:PRINT PTAB(35);"Yes";PTAB(155);"No";
  280. response%=-1 
  281. WHILE response%=-1
  282.  WHILE MOUSE(0)=0:WEND
  283.  mx=MOUSE(1):my=MOUSE(2)
  284.  IF (mx>20 AND mx<80) AND (my>20 AND my<40) THEN
  285.   response%=1
  286.   WHILE MOUSE(0)<>0:WEND
  287.  END IF 
  288.  IF (mx>140 AND mx<190) AND (my>20 AND my<40) THEN
  289.   response%=2
  290.   WHILE MOUSE(0)<>0:WEND
  291.  END IF
  292. WEND
  293. WINDOW CLOSE 4
  294. END SUB
  295.  
  296. SUB introduction STATIC
  297. ' program introduction subroutine
  298. CLS
  299. PRINT
  300. PRINT TAB(4);"J2000 is an interactive AmigaBasic program which can be"
  301. PRINT TAB(4);"used to convert stellar positions, proper motions,"
  302. PRINT TAB(4);"parallax and radial velocity from the standard epoch"
  303. PRINT TAB(4);"B1950 (FK4) to J2000 (FK5). The method used in this"
  304. PRINT TAB(4);"program is based on the algorithm published in the 1985"
  305. PRINT TAB(4);"edition of the Astronomical Almanac, Addendum to Section"
  306. PRINT TAB(4);"B, page X. This technique is an approximation because"
  307. PRINT TAB(4);"it does not account for systematic and individual star"
  308. PRINT TAB(4);"corrections between the FK4 and FK5 epochs."
  309. PRINT
  310. PRINT TAB(4);"J2000 will first request the right ascension, declination,"
  311. PRINT TAB(4);"right ascension proper motion, declination proper motion,"
  312. PRINT TAB(4);"parallax and radial velocity of a stellar object with"
  313. PRINT TAB(4);"respect to the 1950 (FK4) epoch. Each program prompt will"
  314. PRINT TAB(4);"advise you of the proper units to input. If you do not"
  315. PRINT TAB(4);"have a number for proper motion or parallax or radial"
  316. PRINT TAB(4);"velocity, then input '0' for any of these values."
  317. CALL nextpage
  318. END SUB
  319.  
  320. SUB nextpage STATIC
  321. ' select next page subroutine
  322. PRINT
  323. PRINT TAB(13);"< press left mouse button to continue >"
  324. WHILE MOUSE(0)=0:WEND
  325. WHILE MOUSE(0)<>0:WEND
  326. END SUB
  327.  
  328.